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Relaxation of protostellar accretion shocks using the smoothed particle 
hydrodynamics 
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It is believed that protostellar accretion disks to be formed from nearly ballistic infall of the molecular matters in rotating 
core collapse. Collisions of these infalling matters lead to formation of strong supersonic shocks, which if they cool rapidly, 
result in accumulation of that material in a thin structure in the equatorial plane. Here, we investigate the relaxation time 
of the protostellar accretion post-shock gas using the smoothed particle hydrodynamics (SPH). For this purpose, a one- 
dimensional head-on collision of two molecular sheets is considered, and the time evolution of the temperature and density 
of the post-shock region simulated. The results show that in strong supersonic shocks, the temperature of the post-shock 
gas quickly increases proportional to square of the Mach number, and then gradually decreases according to the cooling 
processes. Using a suitable cooling function shows that in appropriate time-scale, the center of the collision, which is at 
the equatorial plane of the core, is converted to a thin dense molecular disk, together with atomic and ionized gases around 
it. This structure for accretion disks may justify the suitable conditions for grain growth and formation of proto-planetary 
entities. 
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1 Introduction 

In the general picture, the gravitational unstable molecu- 
lar cloud cores collapse and form the protostars, accretion 
disks, and outflows. The first model of a spherical symmet- 
ric accretion which flow towards a central object was made 
by Bondi (1952). Years later, Ulrich (1976) modified the 
Bondi's idea by assuming all fluid particles have a certain 
angular momentum with negligible pressure forces at the 
border of the cloud, so that the collapse problem analyzed 
using ballistic trajectories. Since then, many generalizations 
of the core collapse have been made (e.g., Shu 1977, Tere- 
bey, Shu & Cassen 1984, Foster & Chevalier 1993, Galli & 
Shu 1993, Henriksen, Andre & Bontemps 1997, Fatuzzo, 
Adams & Myers 2004, Mendoza, Tejeda & Nagel 2009, 
Bate 2010, Nejad-Asghar2010). As a result of rotating col- 
lapse, the collisions at the equatorial plane, between the in- 
falling matters coming from the northern hemisphere and 
those coming from the southern one, cause the formation of 
strong supersonic shocks. If the shocked gas cools rapidly, 
the result is that material accumulates in a thin structure in 
the equatorial plane, i.e., accretion disk (Hartmann 2009). 

Observations have revealed not only the existence of 
the circumstellar disks around the protostars (e.g., Watson 
et al. 2007, Akeson 2008, Quanz et al. 2010), but also the 
occurrence of jets and outflows (e.g., Hirth, Mundt & Solf 
1997, Wu et al. 2004, Dunham et al. 2010). The outflows 
and jets are ubiquitous together with accretion disks through 
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the collapsing molecular cloud cores (e.g., Furuya, Cesaroni 
& Shinnaga 2011). We cannot directly observe the new- 
born or very young accretion disks and jets because their 
formation places are embedded in a dense infalling enve- 
lope. On the other words, observations cannot directly de- 
termine the real structure and formation history of accre- 
tion disks and outflows. Therefore, theoretical approach and 
simulations are necessary to investigate the formation and 
evolution of them (e.g., Walch et al. 2009, Machida, Inut- 
suka & Matsumoto 2010, Ciardi & Hennebelle 2010, Nejad- 
Asghar 2011). Clearly, consideration of jets and magnetic 
fields can affect on formation and evolution of accretion 
disks (e.g., Hujeirat 1998, Hujeirat et al. 2003, Machida, 
Inutsuka & Matsumoto 2007), but in this study, we turn our 
attention to the phase after formation of of central proto- 
star, and for simplicity neglect the effect of jets and out- 
flows on protostellar accretion shocks. In addition to phys- 
ical evolution of circumstellar disks, the chemical evolu- 
tion from the core to the disk phase is also important (e.g., 
Ceccarelli, Hollenbach & Tielens 1996, Rodgers & Charn- 
ley 2003, Garrod & Herbst 2006, Garrod, Weaver & Herbst 
2008, Visser et al. 2009). Here, we assume that all effects 
of chemical evolution of shocked gas are simplified through 
the appropriate net cooling function. 

Interstellar shocks cover a wide range of parameters: ve- 
locities of 1 — 10 4 km.s~ 1 , pre-shock densities of 10~ 2 — 
10 7 cm~ 3 , and post-shock temperatures of 10 2 — 10 9 K. The 
strength of a shock is indicated by the Mach number which 
can range up to ~ 10 3 , for larger than laboratory shocks 
(McKee & Hollenbach 1980). For an adiabatic strong shock, 
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the jump in density is limited to a factor of four (for a ra- 
tio of specific heats 7 = 5/3), while the temperature of 
post-shock gas can be increased proportional to the square 
of Mach number (e.g., Dyson & Williams 1997). As men- 
tioned above, the high velocity of infalling matters at the 
equatorial plane of a rotating core collapse, leads to strong 
supersonic collisions with great Mach numbers. The high 
temperature of the post-shock gas causes the molecule dis- 
sociation and the atom ionization. Neufeld & Hollenbach 
(1994) examined the physical and chemical processes of 
high-density accretion shocks which associated with the su- 
personic infall of material during the collapse of a molecular 
cloud core to form a protostar. Since the rotation and mag- 
netic fields cause the deflection of infalling matters from the 
protostar so that a shocked disk gas is formed, here we study 
the physical and chemical evolutions of these shocked cir- 
cumstellar disks for understanding the formation of struc- 
ture and substructures through them. 

As a general aspect of star formation, we expect that 
cooling (as a result of chemical and physical changes) of the 
protostellar accretion shocked gas at the equatorial plane of 
collapsing core leads to supply the suitable conditions for 
grain growth and formation of proto-planetary entities. The 
goal of this paper is to investigate this expectation. For this 
purpose, the jump shock (J-shock) structure and the cooling 
time-scale of post-shock gas are presented in section 2. In 
section 3, we use the SPH methodology to investigate the 
time evolution of the strong supersonic shocks. Finally, sec- 
tion 4 is devoted to summary and conclusions. 

2 J-shock structure 

For simplicity, the shock is assumed planar and steady in 
which the deceleration is negligible and there is no ther- 
mal instability in the cooling layer. The jump conditions 
of this shock (J-shock) relate the quantities at an arbitrary 
point behind the shock front to those ahead of it. Conser- 
vation of mass, momentum, and energy across the shock 
front is given by Rankine-Hugoniot conditions (e.g., Dyson 
& Williams 1997) 
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where 7 is the ratio of specific heats, Q is the energy lost 
per unit mass during the shock process, and the equation of 
state is applied as p = (k/ '/j,m#)pT = JCpT. The jump 
conditions (Q~|)-([3]) enable one to solve the J-shock structure. 

We would be interested to consider the collision of two 
gas sheets with velocities vo in the rest frame of laboratory. 
In this reference frame, the post-shock will be at rest and the 
pre-shock velocity is given by v\ — v + v 2 , where v 2 is the 
shock front velocity. Defining the pre-shock sound speed as 



c = \fJC\Ti and Mach number as Mo = vq/c, the equations 
fl)-® lead to 
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where 7" = JC2T2/JC1T1 and TZ = P2I P\ are the relative 
temperature and density of the post-shock gas, respectively. 
Eliminating 7~ between the equations (01 and (0, gives the 
square equation 



AxlZ + A-2 = 0, 
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with coefficients as follows 
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In the general case, the solution of equation (0 can di- 
rectly be found according to four parameters 71, 72, Q and 
Mq; then the relative temperature 7~ is obtained via equa- 
tions (HI or ((5J- The adiabatic shock is a special case with 
Q = 0, which in the strong supersonic collision (Mo — > 
00), the equation © leads to TZ = 4 with assumption of 
72 = 5/3. In this case, the relative temperature is limitless 
as T ~ Mq /3. In fact, this temperature is the maximum al- 
lowed value in a shock process, which is T 2 nax = 3.86^2 x 
10 5 i>q 7 K where v 07 = wo/lOOkm.s -1 (e.g., Hollenbach & 
McKee 1979). The most important parameter in the strong 
supersonic shocks (Mo >> 1) is the energy lost per unit 
mass during the shock process, Q = (n^A/ ' p2mH)tdur, 
where A (erg. cm 3 . s -1 ) is the cooling function at the post- 
shock region with density ri2, and tdur is the duration time 
of the post-shock gas. Accurate determination of the cooling 
time-scale requires specifying the elemental abundance of 
the post-shock region, but a simple estimate can be obtained 
using t coo i s» kT 2 /(n2^)- Eliminating the n^A, we approx- 
imately have Q/c 2 ps T (tdur /tcooi)- If the post-shock gas 
cools rapidly (i.e., t coo \ << tdur), its temperature cannot be 
grater than the molecular dissociation energy, while in slow 
cooling rate (i.e., t coo i >> tdur), the temperature may even 
be grater than 10 4 K causing the atomic ionization process. 

The cooling mechanisms take into account many dif- 
ferent processes that dominate in different ranges of tem- 
perature. We apply the cooling function as outlined in the 
Figure 1 of Heitsch, Hartmann & Burkert (2008) in which 
they used a combination of the rates quoted by Dalgarno and 
McCray (1972) and Wolfire et al. (1995) for T < 10 4 K, 
and the tabulated curves of Sutherland and Dopita (1993) 
for T > 10 4 K. We can express the logarithm of cooling 
function as a piecewise linear function of the temperature 
logarithm as follows 
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Table 1 Parameters for the piecewise linear expression 
of the cooling function, A = Ao(T/Tq)p , with ionization 
degree x% =0.1 and with metallicities corresponding to the 
solar neighborhood. 
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Fig. 1 The cooling time-scale as expressed in a piecewise 
linear form for ionization degree Xi = 0.1 and metallic- 
ities corresponding to the solar neighborhood. The dotted 
curve is the B-spline fit to this piecewise linear expression. 
The temperatures in which the molecules and atoms disso- 
ciate/form are depicted as arrows, and the minimum value 
of cooling time-scale is nominated as trap. 

where A , T and j3 are given in Table 1 for ionization de- 
gree xi =0.1 and metallicities corresponding to the solar 
neighborhood. In this way, the cooling time-scale can be ex- 
pressed in a piecewise form as 

tcooi « — t-T 1 ' 13 , (8) 
n 2 A 

which is shown in Fig.Q~]for n 2 ~ 10 4 cm~ 3 . 

When the strong shocking molecular gases proceed, at 
the first time, the temperature increases near T 2 max so that 
the molecules will be dissociated and the atoms may be ion- 
ized. In the ionized region, there is a trap in which t coo i 



reaches to its minimum value so that the plasma gas cools 
rapidly to recombine the electrons and nucleons. Recom- 
bination becomes important when the post-shock plasma 
cools to 10 4 K; a typical H atomic recombination time-scale 
is approximately t rec sa l/n e a rec « 10 5 /n e (cm -3 )yr, 
where a rec is the total recombination coefficient and n e 
is the number density of electron (Osterbrock & Ferland 
2006). If the process of cooling mode continues, the atomic 
post-shock gas will approach to the molecular gas. The H 2 
molecule cannot form in the gas phase because it is not 
lonely able to radiate away the excess energy of formation. 
The grains, which have been survived in the shock process, 
can thus transfer this extra energy of H 2 formation to their 
phonons (Hollenbach & Salpeter 1971). Although, the de- 
tails of this mechanism are rather complex, a simple esti- 
mation of H 2 formation time-scale is t mo i « l/nsk m « 
10 9 /n 3 ¥ (cm _3 )yr, where k m « 3x 10~ 17 n#(n# + 2uh 2 ) 
cm 3 .s _1 is the hydrogen molecule formation rate on the 
grain surfaces and nn is the hydrogen atomic number den- 
sity (Lequeux 2005). For a typical J-shock structure of a 
protostellar accretion shocks, with n e w nu ~ 10 4 cm -3 , 
we have t rec s» lOyr and t mo x ~ 10~ 3 yr, thus the recom- 
bination occurs slowly while the molecule formation oc- 
curs very fast. Although, we find some insights for J-shock 
structure and its relaxation process, but the shock process 
is rather a complex non-equilibrium thermal and chemical 
mechanism. Thus, in the next section, we use the SPH sim- 
ulation to study the time evolution of processes in the strong 
shocking molecular gases. 

3 Investigation by SPH 

We consider the head-on collisions of two molecular gas 
sheets. An initial negative velocity is given to the particles 
with a positive ^-coordinate, and a velocity in the oppo- 
site direction to those with a negative x-coordinate. Com- 
positions of the sheets are assumed as global neutral which 
consists of a mixture of atomic and molecular hydrogen 
(X sa 0.75), helium (Y sa 0.25), and traces of CO and 
other rare molecules. The mean molecular weight is ini- 
tially 1/fj, = X/2 + y/4 sa 0.4375 for molecular case 
(T < 10 2 K), and in the post-shock region for simplicity 
is assumed to be 1/p, = X + Y/4 s» 0.8125 for atomic case 
(10 2 K < T < 10 4 K), and 1/fj, = 2X + 3Y/4 sa 1.6875 
for ionized case (T > 10 4 K). The same procedure is fol- 
lowed for the ratio of specific heats: 7 = 7/5 for diatomic 
molecular gas and 7 = 5/3 for atomic and ionized cases. 
The chosen physical scales for length and time are [I] = 
3.0 x 10 14 cm = 20AU and [t] = 3.0 x 10 7 s = lyr, respec- 
tively, so the velocity unit is approximately 100 km.s -1 . 
The gravitational constant is set G = 10~ 12 [Z] 3 . [t] — 2 . [ttt,] — 1 
for which the calculated mass unit is [to] = 4.5 x 10 23 g. 
Consequently, the physical scales for density and energy are 
[p] = 1.6 x 10- 20 g.cm- 3 = 10 4 H.cm- 3 and [e] = 4.5 x 
10 30 erg, respectively. Two equal one dimensional molecu- 
lar sheets with extension x = 0.1 [Z] is considered, which 
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have initial uniform density and temperature of 10 4 cm 3 
and ~ 10K, respectively. 

The SPH method is well suited to address unbound as- 
trophysical problems, especially the behavior of gases sub- 
jected to compression (Rosswog 2009). A worthy review 
of the SPH methodology and its applications can be found 
in Monaghan (2005). In this method, fluid is represented by 
N discrete but extended/smoothed particles (i.e. Lagrangian 
sample points). The particles are overlapping, so that all 
the involved physical quantities can be treated as continu- 
ous functions both in space and time. Overlapping is rep- 
resented by the kernel function, W a b = W(r a — r^, h a b), 
where h a b = (h a + h},)/2 is the mean smoothing length 
of two particles a and b. The density is estimated via usual 
summation over neighboring particles, 



Pa 



b 



m b W ab , 



(9) 



the acceleration equation in the one-dimensional usual sym- 
metric form is 



d-Vg 

dt 



u Pa 



Pa ^Pb_ 
Pi 



n ab ) 



dWab 
dx a 



(10) 



and the SPH equivalent of the energy equation is 



dUg 

~~dt~ 
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where u a — - 1 _ 1 lC a T a is the thermal energy per unit mass, 
and H a b is the artificial viscosity between particles a and b 



0. 



if Vab*ab < 0, 

otherwise, 



(12) 



where \ a b = \ a — \b and r a b = r a — Tb are relative velocity 
and the distance of particles, p a b = \ (p a + p b ) is an average 
density, v S i g = (c a + Cb)/2 is signal velocity where c a and 
c b are the sound speed of particles, and /i* fc is defined as its 
usual form 

1 



Vab 



(13) 



hab r 2 Jhl b + T? 2 
with ij ~ 0.1. 

To prevent unphysical solutions with inter-particle pene- 
tration and unwanted heating, we use a* and j3* in the form 
of variables with respect to time, 

</n " - al ~ amin +z a S a , (14) 



and 



dt 

dPl 
dt 



(15) 



respectively, where S a = max(-|^, ^ - ^ - ^f, 0) is 
the restricted source term, r a = h a /(Cv S i g ) with 0.1 < 
C < 0.2 is the decay time-scale, and the parameters z a 
and are chosen to regulate the effect of source term so 




(a) 



0.00 

x 



(b) 

Fig. 2 (a) Temperature and (b) density of the adiabatic 
shock (Ao = 0) in the head-on collision of two sheets with 
initial Mach number Mq = 500. The solid curves are de- 
rived with variable viscosity (fT4l and ([TBI , while the dot- 
ted ones are from the common artificial viscosity of Mon- 
aghan (1989) with a* = 1 and (3* = 2. The analytical re- 
sults of strong supersonic adiabatic shocks are depicted by 
arrows. 



that the heat production and post-shock oscillations are con- 
trolled in the numerical simulations (Nejad-Asghar, Khesali 
& Soltani 2008). Here, we choose a m i n — 1, f} m in = 2, 
C = 0.2, z a — 3, and zp — 0.01 for the best smoothed re- 
sults of simulations (see, e.g., Fig. [2j. Clearly, comparison 
between numerical and analytical results in Fig. |2] verifies 
the accuracy and physical consistency of the employed nu- 
merical tool. 

The width of the post-shock region increases by time, 
until it reaches to all extensions of the simulated sheets, or 
the temperature of center of the shocked region cools less 
than 10K. For computer experiments which are considered 
here, with initial extension x = 0.1 [Z], we stop the simula- 
tion when 90% of total SPH particles enter to the post-shock 
region or the temperature of center of the shocked region 
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Fig. 3 Post-shock temperature at the simulation epoch for 
different initial Mach numbers. The adiabatic case, T2 oc 
Mq , is shown by dash-line. 

cools less than 10K. In this simulation epoch, the cooling 
rate affects on the post-shock temperature as shown in Fig. [3] 
for various Mach numbers. In this figure, the relationship 
T2 = j^-^Mq, as mentioned for the strong supersonic adi- 
abatic shocks, is depicted by dash-line. We see the cooling 
rate causes to decrease the value of post-shock temperature. 
The effect of cooling rate appears more clear in the Mach 
numbers which cause to set the post-shock temperature in 
the trap region as outlined by Fig. [2] Since for larger Mach 
numbers, the simulation epoch (i.e, the time in which 90% 
of SPH particles of our simulation enter to the post-shock 
region or the temperature of center of the shocked region 
cools less than 10K) is shorter, the post-shock temperature 
reaches asymptotically to the adiabatic case. 

In the general case, the high temperature of the post- 
shock region lead to splash it into the medium. But, in the 
accretion shock, the infall matter confides the shocked gas 
so that it may cool and end up with an accretion rotating 
disk. For finding the relaxation time, we use the equation 
([TlT l with v a b ~ 0, thus, we have 



1 „ dT a 

A-a — ; — 

1 dt 



A 



7a - 1 at (MaTOff) 

The temperature of the accreted shocked gas can be ob- 
tained by integrating the equation ([Tol l. The result is in a 
piecewise change from T a \ (at t{) to T a2 (at t 2 ) as follows: 

^ (l-/3)( 7a -l) A 



T a2 = {T^- 



tc a 



Pa_ 

-L n 



(/i Q ™ff) 2 

(ta-ti)}^, 



(17) 



which is shown in Fig.@] with assumption of p a s» 4. 

4 Summary and conclusions 

Molecular cloud cores are rotating so that in the process 
of collapse and protostellar formation, the infalling matters, 



Fig. 4 Relaxation of the temperature of the post-shock gas 
in the head-on collision of two sheets. 



which arrive at the equator, collide and dissipate the kinetic 
energy of motion perpendicular to the equatorial plane so 
that an accretion disk can be formed. The speed of infalling 
matters at the equator is so high that the strong supersonic 
shocks appear, and the temperature of post-shock gas is in- 
creased so causing to dissociate the molecules and ionize 
the atoms. In adiabatic strong supersonic shocks, the density 
of post-shock region is about fourfold of initial density, and 
the temperature is increased proportional to the square of 
the Mach number. On the other hand, the cooling processes 
of the post-shock gas can decrease the temperature so that 
the ionized gases can be recombined to form the atoms and 
molecules, if the duration time-scale of the post-shock re- 
gion is longer than the cooling time-scale. The suitable cool- 
ing function for the post-shock gas (Table 1 ) showed that the 
cooling time-scale has a minimum which is at the ranges of 
ionizing temperature (Fig.[T). This minimum, which is nom- 
inated as a trap, causes to cool the post-shock gas in a fast 
rate. Thus, if the initial speed of colliding matters is so high 
that the temperature of post-shock gas settles in this trap 
region, it will quickly cool and electrons recombine with 
nucleons. 

For investigating the time evolution of the post-shock 
gas in the strong supersonic collisions, which occurs at pro- 
tostellar accretion shocks, we used the SPH method with 
the variable artificial viscosity to obtain the more smoothed 
results of simulations (Fig. The simulations of strong 
shocks with great Mach numbers show that the temperature 
of the post-shock gas is quickly increased near to T™ ax , 
which is for the adiabatic case, and gradually decreased by 
the time according to the cooling rate. Decrease of tempera- 
ture of the post-shock gas at the simulation epoch is shown 
in Fig|3]for different values of the Mach numbers. We see 
from this figure that the decreasing rate of temperature is 
very fast in the trap region as depicted by Fig. Q] 

The simulations show that the temperature of center of 
the post-shock gradually decreases, while in the ridges, it 
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stays about T™ ax because of continuously infall of matters. 
Temperature decrease of the central region leads to increase 
of its density in an isobaric manner so that an accretion thin 
disk can be formed. Thus, over the time, center of the col- 
lisional infalling matters (i.e., equatorial plane) converts to 
a dense molecular thin disk with an atomic envelope and 
ionized gas which comes from strong shocks of continuous 
infalling matters. The time in which this structure occurs, 
depends on the Mach number that is shown in Fig. [4] The 
temperature range of the trap, which leads to fast cooling 
of the post-shock gas, is clearly seen in Fig. |4] too. Thus, 
we see that the cooling processes of the post-shock gas in 
the protostellar accretion shocks can lead to the formation 
of an accretion thin molecular disk at the equatorial plane, 
in a convenient time-scale. This dense molecular thin disk 
is appropriate for grain coagulation and formation of proto- 
planetary entities. 
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